function L=CENP_E_tracking(MT,T,savefolder) %T -is the upper time limit for calculation

%initial conditions------------------------------------------------------
XMTend=1000; %initial MT tip coordinate, nm
Xhead1=8; % initial motor domain coordinate, nm
Xhead2=16; % initial motor domain coordinate, nm
Xtail=0; % initial tail domain coordinate, nm
head_detached=0;
tail_detached=0;
      
%model parameters--------------------------------------------------------------
if ~((MT=='d')||(MT=='p'))
      disp('Warning: MT polymerization status was not specified correctly, so it was set to "depolymerization" by default');
      MT='d';
end
if MT=='p'
    p_on_tubulin=20.5; %sec-1
    p_off_tubulin=17.4; %sec-1
elseif MT=='d'
    p_on_tubulin=0;  
    p_off_tubulin=28.8; %sec-1
end
tail_diffusion=1.6*10^6;
v_motor=333.333;%nm/sec
p_head_detached=0.048; % sec-1
p_tail_detached=2; % sec-1
p_reattached=2600; % sec-1
pinhib=0.0625; % sec-1

%technical parameters------------------------------------------------
dt=0.00002; %sec, time resolution of the calculation. 20 usec is recommended!
dt_for_plots=0.01; %sec, time resolution of the output graphs

%--------------------------------------------------------------------------
left=1; % variable used to describe hand-over-hand stepping
right=0; % variable used to describe hand-over-hand stepping
tip_tracking_started=0; %flag for determining the moment of tip-tracking onset
i=1;
it=1;
N=floor(T/dt);
TAIL=zeros(N,1);
HEAD1=zeros(N,1);
HEAD2=zeros(N,1);
TIP=zeros(N,1);
time=zeros(N,1);

% verifying that initial conditions are set correctly -------------------------------
if abs(Xhead1-Xhead2)>8
      disp('Warning: Wrong distance between the heads. Heads positions set to default');
      Xhead1=8; 
      Xhead2=16; 
end
if (abs((Xhead1+Xhead2)/2-Xtail)>230)||(Xtail>Xhead1)||(Xtail>Xhead2)
      disp('Warning: Wrong distance between the tail and the heads. Positions set to default');
      Xtail=(Xhead1+Xhead2)/2-16; 
end
if (XMTend<=Xhead1)||(XMTend<=Xhead2)||(XMTend<=Xtail)
      disp('Warning: Wrong position of CENP-E relative to MT end. Positions set to default');
      XMTend=max([Xhead1 Xhead2 Xtail])+1000; 
end
%--------------------------------------------------------------------------------------------

while i<N
    if i==floor(it*N/10)
          disp(['progress: ',num2str(it*10),'%']); %shows progress relative to maximal time of calculation
          it=it+1;
    end
    %MT dynamics
    randnum=rand(1);
    tub_falloff=0;
    if randnum<1-exp(-p_on_tubulin*dt)
         XMTend=XMTend+8;
    end
    if randnum<1-exp(-p_off_tubulin*dt)
         tub_falloff=1;
         XMTend=XMTend-8;
    end
    
    %Head and tail reattachments         
    randnum=rand(1);
    if randnum<1-exp(-p_reattached*dt)
       if (head_detached==1)&&(tail_detached==0)
             Xhead1=Xtail+38*randn(1)+135; 
             %38 and 135 above are STD and mean of head-to-tail
             %distribution, obtained by fitting WLC model with persistence length=26nm
             
             Xhead1=(round(round(Xhead1)/8))*8;%making sure that attachment occurs to the binding site
             while Xhead1>XMTend; %sampling until the head's coordinate becomes within MT length
                   Xhead1=Xtail+38*randn(1)+135;
                   Xhead1=(round(round(Xhead1)/8))*8;%making sure that attachment occurs to the binding site
             end
             right=1;
             Xhead2=Xhead1-8;
       end
       if (head_detached==0)&&(tail_detached==1)
           Xtail=Xhead1+38*randn(1)-135;
           while Xtail>XMTend;
                 Xtail=Xhead1+38*randn(1)-135;
           end
       end
       head_detached=0;
       tail_detached=0;
    end
    
    %Heads ATP-dependent motion
    randnum=rand(1);
    if head_detached==0
       if (XMTend-Xhead1<=8)||(XMTend-Xhead2<=8)
          if tip_tracking_started<1
               tip_tracking_started=1;
               Xstart=XMTend;
          end
               
          left=0;
          right=0;
       else
           if Xhead1>Xhead2
              right=1;
              left=0;
           else
              right=0;
              left=1;
           end
       end
       if randnum<1-exp(-v_motor/8*dt)
            Xhead1=Xhead1+left*16;
            Xhead2=Xhead2+right*16;
            temp=left;
            left=right;
            right=temp;
       end
    end
    
    %head's falloff with detaching tubulins
    if left+right<1%(XMTend-Xhead1<8)||(XMTend-Xhead2<8)
        if tub_falloff==1
            head_detached=1;
            if tail_detached<1
                    Xhead1=Xtail;
                    Xhead2=Xtail;
            end
        end
    end
    
    %head falloff by itself
    randnum=rand(1);
    if randnum<1-exp(-dt*p_head_detached)
        head_detached=1;
        if tail_detached<1
                Xhead1=Xtail;
                Xhead2=Xtail;
        end
    end
    
    %tail behavior
    if head_detached<1
        Xtail=min([Xhead1, Xhead2])+38*randn(1)-135; %tail stays near the head with distance determined from WLC model
    else
        dXtail=randn(1)*sqrt(tail_diffusion*2*dt);
        if Xtail+dXtail<XMTend
           Xtail=Xtail+dXtail; % free 1-d diffusion along MT
        else
           Xtail=2*XMTend-Xtail-dXtail; % tail reflection from MT end
        end
    end
    
    %tail detachment
    randnum=rand(1);
    if randnum<1-exp(-dt*p_tail_detached)
        tail_detached=1;
        Xtail=Xhead1;
    end
    
    %tail detachment with falling off tubulin
    if (XMTend-Xtail<8)&&(tub_falloff==1)
            tail_detached=1;
    end

    %self-inhibition    
    if tail_detached+head_detached<1
        randnum=rand(1);
        if randnum<1-exp(-pinhib*dt)
           tail_detached=1;
           head_detached=1;
        end
    end
    
    TAIL(i)=Xtail;
    HEAD1(i)=Xhead1;
    HEAD2(i)=Xhead2;
    TIP(i)=XMTend;
    time(i)=dt*i;
    if tail_detached+head_detached==2 %complete detachment
          break;
    end
    i=i+1;
end
i_end=i;
num=1:i_end;
TIP=TIP(num);
TAIL=TAIL(num);
HEAD1=HEAD1(num);
HEAD2=HEAD2(num);
time=time(num);

%calculate CENP-E's tip-tracking run length
if tip_tracking_started>0
      L=round(abs(0.5*HEAD1(end)+0.5*HEAD2(end)-Xstart));
else
      L=0;
end
disp(['tip tracking run length = ',num2str(L),' nm']);


if dt_for_plots<dt %check if the default plot resolution can be obtained in this calculation
    disp('Warning: time step for calcs is set to be lower than default time step for graphs. So the time step for graphs will be increased.');
    dt_for_plots=dt; 
end
everyN=floor(dt_for_plots/dt);
TIPdown=downsample(TIP,everyN);
TAILdown=downsample(TAIL,everyN);
HEAD1down=downsample(HEAD1,everyN);
HEAD2down=downsample(HEAD2,everyN);
TIMEdown=downsample(time,everyN);

%plot everything in a way similar to a kymograph
figure 
plot(TIPdown,TIMEdown,'r.-',TAILdown,TIMEdown,'g.-',HEAD1down,TIMEdown,'k.-');
xlabel('Coordinate along MT axis, nm');
ylabel('Elapsed time, sec');
legend('MT tip','TAIL','Head1','Location','EastOutside')
hold on
axis ij

%save_data
dlmwrite([savefolder,'\time_tip_tail_head1_head2.txt'],[TIMEdown TIPdown TAILdown HEAD1down HEAD2down],'delimiter','\t','newline','pc');


